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ABSTRACT 

We present investigations of the interaction between the relativistic, precessing jets of the 
microquasar SS 433 with the surrounding, expanding Supernova Remnant (SNR) shell W 50, 
and the consequent evolution in the inhomogeneous Interstellar Medium (ISM). We model 
their evolution using the hydrodynamic FLASH code, which uses adaptive mesh refinement. 
We show that the peculiar morphology of the entire nebula can be reproduced to a good 
approximation, due to the combined eff'ects of: (i) the evolution of the SNR shell from the 
free-expansion phase through the Sedov blast wave in an exponential density profile from 
the Milky Way disc, and (ii) the subsequent interaction of the relativistic, precessing jets of 
SS 433. Our simulations reveal: (1) Independent measurement of the Galaxy scale-height and 
density local to SS433 (as hq = 0.2cm"^,Zd = 40 pc), with this scale-height being in excel- 
lent agreement with the work of Dehnen & Binney. (2) A new mechanism for hydrodynamic 
refocusing of conical jets. (3) The current jet precession characteristics do not simply extrap- 
olate back to produce the lobes of W 50 but a history of episodic jet activity having at least 3 
diff'erent outbursts with diff'erent precession characteristics would be sufficient to produce the 
W 50 nebula. A history of intermittent episodes of jet activity from SS 433 is also suggested in 
a kinematic study of W 50 detailed in a companion paper (Goodall et al, MNRAS submitted). 
(4) An estimate of the age of W50, and equivalently the age of SS 433 's black hole created 
during the supernova explosion, in the range of 17,000 - 21,000 years. 

Key words: Hydrodynamics - supernovae: general - ISM: jets and outflows - ISM: kine- 
matics and dynamics - ISM: supernova remnants - X-rays: binaries 



1 INTRODUCTION TO THE SS 433-W 50 COMPLEX 

The paradigmatic microquasar SS 433 has become widely known 
since its discovery as the first stellar source of relativistic jets 
within our G alaxy ( [Margon et aL| 1 979 1 [Fabian & Rees|1979|[Mir] 
|grom|l97"9l l. After three decades of researct[^ the system has been 
widely identified as a high mass X-ray binary (HMXB) system. The 
masses of the constituent bodies have been estimated using a vari- 
ety of techniques (Crampton & Hutchings 1981 , D'Odorico et al.| 



around SS 433 has been observed via optical spectros copy (|Blun 



1991 , Collins & Scher 2002, Gies et al. 2002, Hillwig et a l.|2004[ 
Fuchs et al. 2006, Lopez et al. 2006 , Blundell et al. 2008 ) yielding 
some varied results. In general the observations seem to favour the 
scenario of a black hole in orbit around a more massive companion 
star ( |Fabrika|2004j . Recently, the presence of a circumbinary ring 
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dell et al.|2008| ) and also near infra-red spectroscopy ( [Perez M. & 



Blundell[20 09X and these observations determine the mass internal 



to the ring to be 40Mo, of which 16Mo is attributable to the com 
pact object and its accretion disc. Thus, it is believed that SS 433 
hosts a black hole (BH) rather than a neutron star. The BH's ac- 
cretion disc is presumably fed by gas from a strong stellar wind 
or Roche-lobe-like overflow from the companion s tar with a rate 
of 7 X 10"^ < Mtransfer < 4 X 10"^ Mq ycar-^ ( Ki ng et al.|2000[ ), ac- 
companied by two oppositely directed relativistic jets wh ich are 
thought to eject material at a rate Mjet > 10~^ M© year~^ ( [Begel 



[man et alT I 1980) into the surrounding interstellar medium. The 
disc -jet system exhibits precessional motion, which is roughly de- 
scribed by the kinematic model of SS 433 (Milgrom,,1979, ,Abell[ 
[& Margon|1979| and is readily observed via the periodic Doppler 
shifting of optical emission lines from the jets. By fitting the kine- 
matic model to spectroscopic data spanning 20 years, |E ikenberry| 
[et al.| ( |2001 ) determined the following best-fit kinematic parame- 
ters: fpi-ec = (162.375 + 0.011) days for the precession period, ySjet 
= (0.2647 + 0.0008) for the jet speed, ^prec = (20°. 92 + 0°.08) for 
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The location of W 50 within the Milky Way 
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Figur e 1. This figure describes the geometry of the SS 433-W 50 system, (a) This mosaic was created using archival data from the GBT6 survey ( [Griffith et al.] 
|l995) , showing the location of W 50 with respect to the Milky Way Galaxy plane. The vector GD points normally towards the Galactic disc from SS 433, and 
^GD = 19.7° is the angle this vector subtends with the mean jet axis. The plate scale in (a) applies to objects at the same distance as SS 433 of 5.5kpc (1 arcmin 
^1.6 pc). (b) The famous VLA mosaic of W50 provided by Dubner et al. (1998) has been transformed here such that the x-axis is coincident with the mean 
axis of SS 433's jet precession cone. The solid green lines indicate the current jet precession cone in SS 433 which subtends an angle of Op^ec with the mean 
jet axis. SS433's western jet-cone points in the direction +x, and the eastern jet-cone in -x. In this frame the Right- Ascension ordinate a makes an angle of 
172° with the x-axis (6a = 8°). An example of the Galactic density profile used in our model with Z^=40pc is indicated by the dashed black lines. Adjacent 
dashed black lines are separated by (Zd In 2) parsec corresponding to changes by a factor of 2 in density, and the density profile is normalised to «o=l particle 
cm~^ at SS433. Another microquasar, GRS 1915-1-105 is also present in the field of view of (a) and its coordinates are indicated. However, GRS 1915-1-105 
does not seem to be encapsulated within a well defined large-scale nebulous structure in the way that SS 433 is with W 50. 



the precession cone half-angle, and Oinc = (78°. 05 + 0°.05) for the 
inclination of the mean jet axis to our line of sight, although [Blun-| 
I dell & Bowler ( |2005|) found that these properties vary with time. 
'Goranskii et al. ( 1998 ) conducted a multi-colour photometric study 
of SS 433 over a similar time period of 20 years, and used Fourier 
analysis to extract the orbital period of Porb = 13.082 days, and a 
nutation period of Pnut = 6.28 days. 

In the radio regime, more than two complete precessional cy- 
cles of the jets have been spatially resolved ( [Blundell & Bowler] 
[2004| l, providing confirmation of the jet speed as ~0.26c and en- 
abling an independent calculation of the distance to SS433 as 
<^ss433 = 5.5 + 0.2 kpc. Using Very Long Baseline Interferometry 
(VLB I), the jet has been resolved at milliarcsec scales into individ- 
ual fireballs or bolides of ejecta ( [Vermeulen et al.|1987l|Paragi et al.| 
[1999bj ), and the presence of an equatorial ''ruff" outflow from the 
accretion disc has become apparent ( [Paragi et al.|1999a| |Blundell| 
'et al. 2001 ). Furthermore, recent optical spectroscopic studies have 
revealed that these discrete fireballs of jet ejecta are optically thin, 
and expanding in their own rest frame at approximately 1 % of the 
jet speed ( [Blundell et al.|2007| ). 

SS 433 is embedded within the W 50 nebula, as shown in the 
radio mosaic of [Dubner et al.[ ( p^98] ). The remarkable phenomenol- 
ogy of W 50 (Figure[T]3) spans 208 pc along its major axis, featuring 
a circular region which is believed to be the remains of an expand- 
ing supernova remnant (SNR) shell, with two elongated lobes to the 
east and west of the circular shell. The circular shell is centred upon 
the position of SS 433 with an off'set of just 5 arcminutes ( [Lockman] 
[et al.|2007] ). The distance to W50 is difficult to constrain using the 
kinematics of the HI gas within the nebula because of its extent 
and inhomogeneity, and due to the high level of turbulent confu- 
sion along the line of sight near the Galactic plane. However, recent 



radio observations of HI in absorption and emission now confirm 
that the distance to W 50 is consistent with 5.5 kpc ( Lockman et al.[ 
[2007] ), which places it at the locality of SS 433. Hence we consider 
the details of the local environments of SS 433 and W 50, which 
are the same distance from us, such that their mutual evolutions are 
connected. 



2 PREVALENT QUESTIONS PERTAINING TO THE 
SS 433-W 50 COMPLEX 



The detail revealed by the famous radio mosaic of W 50 ( [Dubner[ 
[et al. |1998 ) presented a number of puzzling questions concerning 
its formation. We endeavour to address these issues through our 
models and will later discuss each one (see ^6] iTland ^8). 



2.1 East-west lobe asymmetry and annular structure 

The lobes share an axis of symmetry with SS 433's mean jet axis, 
and exhibit a fascinating annular structure (most noticeable in the 
eastern lobe, Figure[T]3). This has led to the hypothesis that SS 433's 
jets are responsible for the formation of W50's elongated lobes. 
There is a noticeable asymmetry in the extents of the east and west 
lobes of W 50 (Figure[T]3), and this can be characterised by the ratio 

-^East 



Pratin — 



1.40. 



(1) 



Due to the inclination of the jet- axis with respect to our line of 
sight, a small geometrical correction sin Oi^c ~ 1 is implicit in the 
measured values Pnast and Pwest (but not in the ratio Pratio). where 
^inc ^ 80° ( .Hjellming & Johnston J 98T|[Eikenberry et al.[2001| ). As 
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a result the physical extent of each lobe is larger than measured, by 
a factor 1/sin ^inc- This asymmetry between the eastern and western 
lobes is thought to be due to the exponential ISM density gradi- 
ent towards the Milky Way disc, and the validity of this theory is 
explored in ^7.1| of this paper. 

2.2 SS 433's latency period 

We can calculate an in vacuo jet travel time required for the extent 
of the eastern lobe to reach /?East based on the jet speed: 



^East . cr\c\ 

^vacuo - - — : - 1500 years 



(2) 



where ^vacuo is only a lower-limit estimate for the actual jet travel 
time; in reality the travel- time would be longer due to retardation 
of the jets in the presence of the ISM. This estimate is noticeably 
smaller than the time taken for a typical SNR to reach a radius of 45 
parsecs (the radius of the circular region in Figure [T]3), suggesting 
that SS 433's jets were initiated relatively recently in the history 
of W 50. We refer to the period after the supernova explosion of 
SS 433's progenitor and before the ignition of jets, as SS 433's "la- 
tency" period. The reason for this latency period is currently not 
well understood, and we endeavour to estimate the duration of this 
period in ^in order to constrain the possible scenarios. 

2.3 SS 433's jet persistence and precession-cone angle 

SS 433 is often noted as being somewhat unique among the micro- 
quasar class due to the apparent persistence of its radio jets. Other 
microquasars typically exhibit a more intermittent behaviour, tem- 
porarily increasing in luminosity with each episodic jet outburst, 
before returning to a state of relative quiescence. These jet ejec- 
tion episodes in microquasars characteristically last for one or more 
weeks, and recur every one or more years. Successive jet outbursts 
can be observed to happen with substantially different properties, 
for example in 1997 a flare from Cygnus X-3 was observed with 
jet speed V1997 > 0.81c and precession cone angle 0^1997 < 12° 
( [Mioduszewski et al.|2001| ), while a subsequent flare was observed 
in 2001 with jet speed V2001 = 0.63c and precession cone angle 
4^1001 - 2.4° ( [Miller- Jones et al.|2004| ). We must consider the pos- 
sibility that SS 433's jets are not constantly active, but rather that 
SS 433 also exhibits jet outburst episodes, albeit of a much longer 
duratiorj^than has been observed in other known microquasars. 

2.4 Previous simulation work pertaining to W 50 

The range in physical sizes involved in the SS433-W50 system 
makes it an extremely difficult object to model numerically, how- 
ever several authors have confronted the challenge of simulat- 
ing the jet-SNR interaction (notably Kochanek & Hawleyl ( |1990| ; 
[Velazquez & Raga| { [2000| ; [Zavala et al.|(|2008| )). Although a more 
detailed comparison will be given in |9.3[ we briefly introduce 
some of the early approaches here. 

An early hydrodynamical study of SS 433's jets was carried 
out by Kochanek & Hawley[ ( [1990 ), who considered various forms 
of hollow cylindrical and conical jets, as well as filled conical jet 
models. The effect of "self interaction" of the jets with previous 
ejecta and with the jet cocoon itself is mentioned in their paper, 

^ In other words, we are currently witnessing (and have been since its dis- 
covery) SS 433's most recent outburst of jet activity. 



and in several of the references therein. They use an axisymmetric 
model, and report that their hollow conical jets are neither recol- 
limated, nor refocused. They discuss how their results are in dis- 
agreement with a model for jet recollimation proposed by [Eichler| 
|1983^, in which precessing jets were suggested to experience a fo- 
cusing effect caused by the pressure from the ambient ISM through 
which they propagate. 

[Velazquez & Raga[ pOOO] ) recognised the physical attractive- 
ness of the jet-SNR interaction scenario. They modelled the W 50 
nebula by instantaneously imposing an analytical Sedov profile 
upon the ISM to mimic the supernova explosion, upon which the 
relativistic jets of SS 433 are later incident. Due to the computa- 
tional constraints at the time, they modelled only one quadrant of 
the SNR in two dimensions, with the simplification of a scaled- 
down version of W50 of radius lOpc and a uniform density back- 
ground ISM. The system evolution was closely followed until 2700 
years after the jets were switched on, such that the simulated jet- 
SNR lobe extends to approximately 33pc from the SNR origin. 
They follow the example of [Clarke et al.[(T989| ) who describe a 
neat model for computing the synchrotron emission from 2-D sim- 
ulations, by assuming a toroidal magnetic field configuration about 
the jet axis, and revolving the resulting emission about In radi- 
ans to create the projected 3-D emission. The hydrodynamics of 
the jet-SNR was further studied in the 3-d simulations of [Zavala] 
[et al.| (2008 ), who invoke both hollow conical jets to simulate the 
precession in SS 433. They also conduct a small number of tests to 
investigate the physical parameters such as jet mass-loss rate and 
the radius of the SNR at the time when the jets switch on. The pos- 
sibility of a time- variable jet precession angle is also tested in their 
investigations, and their results show a qualitatively similar mor- 
phology to that observed in W 50, in the sense that a spherical SNR 
shell is present with lateral extensions in the form of jet lobes, of 
arbitrary size. 

Despite the successes of these pioneering works, there are sev- 
eral invalid assumptions made in the models which have significant 
effects upon the results, and can be improved upon. The physical 
significance of the morphology present in W50 has (so far) not 
been adequately addressed with relevance to the pressing questions 
(see The inconsistencies in these previous works will be dis- 
cussed more fully in section g9.3| and we endeavour to show that it 
is possible to produce a model of the SS 433-W50 interaction that 
is consistent with observational data, without imposing any unrea- 
sonable assumptions. 

We emphasise the importance of including relevant observa- 
tional data in simulation models, such as the actual sizes of the SNR 
and jet lobes, if we are to provide useful constraints on the system. 
We incorporate all of these observational constraints (with the ex- 
ception of orbital motion) into our hydrodynamical model, to give 
an accurate representation of the dynamics of SS 433's jets. Using 
the observed parameters as a basis for our models of the SS 433- 
W 50 complex, we consider three possible evolutionary scenarios 
in detail in order to see if we can explain the origins of this ex- 
tremely large and rare nebula. The results of these three scenarios 
are discussed in detail in |6] ^and |8] but their mechanisms can 
be summarised as follows: (i) a supernova blast wave expanding in 
the presence of a realistic exponential baryon density background 
appropriate to SS 433's location in the Milky Way (ii) the interac- 
tion of SS 433's jets with the ISM only (no supernova), and finally 
(iii) the interaction of SS 433's jets with an evolved SNR shell from 
(i). 

However, our intention is not solely to replicate the morphol- 
ogy of SS433-W50 complex, but to explore the effects of impor- 
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tant astrophysical parameters such as the jet mass-loss rate, super- 
nova blast energy, and the ambient ISM density, which are among 
input parameters to our supernova and jet models. The distinctive- 
ness of the W 50 nebula produced through these interactions, allows 
us to place useful constraints upon these system parameters. 



3 OUR NUMERICAL SIMULATIONS 

To conduct our simu lations, we utilise t he well-known hydrody- 
namic FLASlfjcode (| Fryxell et al. 2000| . The code implements an 
adaptive mesh refinement (AMR) method, a higher order Godunov 



method known as the Piecewise Parabolic Method (PPM) ( Colella 
|& Woodward] 19 84 j . FLASH is noted for its flexibility, and several 
case- specific enhancements were made by the authors; most no- 
tably the implementation of an exponential Galactic density back- 
ground and our supernova explosion and jet model modules, which 
are described in this section. 

High resolution VLBI imaging of the central engine in SS 433 
reveals jet structure on scales of a few tens of AU, a factor of 10^ or 
so smaller than the extent of W 5(j^ The need for a high resolution 
representation of the core region of SS 433, and the large-scale field 
of view (FOV) necessary to enclose W 50, makes for a challenging 
computational problem. 

Figure[T]shows the location and orientation of W 50 within the 
Milky Way and defines the axes of the simulation domain; the x and 
y axes lie in the plane of the page, and the z axis points out of the 
page, with SS 433 at the coordinate origin. To maximise our reso- 
lution in the regions of interest, we simulate a thin slice through the 
centre of the SS 433-W 50 complex, with grid dimensions (Lx X Ly) 
= (230 X 115) parsecs. The grid is initially split into (/ixb x fZyb) 
blocks before any "refinement" of the grid takes place, and it is 
these blocks that are subject to refinement. In order to ensure that 
all elements on the grid are square rather than rectangular, we set 
«xb=16 and f7yb=8. Our chosen domain size ensures that neither 
jets nor the SNR reach the domain boundaries, such that outflow 
boundary conditions can be suitably applied. Our limiting resolu- 
tion depends on the maximum level of refinement (//max) allowed, 
such that: 



dx - 



[Wxb X 2''max-l 



= dy = 



Ly 



[Wyb X 2'7max-l ] 



(3) 



The domain thickness in the z direction is not constant; each 
cell across the simulation domain has a thickness in z equal to its 
length in x (which depends on the cell refinement level), and so the 
domain is eff'ectively 2D (once cell thick in z). 

We set the maximum refinement level to //max = 10, resulting 
in a maximum resolution of 6x-6y- 0.028 pc = 8.66x10^^ metres 
on our grid, or ^ 1 arcsecond on the sky for objects at a distance 
of Jw5o = 5.5kpc. This resolution corresponds to about three times 
the synthesised beam- width from the famous VLA observation of 
SS433's corkscrew jets ( [Blundell & Bowler|2004| ). For simplicity 
we refer to the maximally refined cells as pixels (cells for which 



^ceii = ^max), and duc to the nature of AMR, every cell in the sim- 
ulation domain therefore has sides equal to an integer power of 2 
times the pixel size. Where present, radiative cooling refers to the 
standard cooling functions provided with FLASH2.5. The standard 
radiative cooling in FLASH2.5 assumes an optically thin plasma, 
where the cooling function h{T) is based upon models of the en- 
ergy losses from the transition region of the solar corona |Rosner| 
|et al.|1978 ) and the chromosphere (Peres et al. 1982 ), and the re- 
sulting cooling function is given as a piecewise-power law approx- 
imation (refer to the user guide available on the FLASH website). 
We acknowledge that the cooling function adopted by FLASH2.5 
assumes that the ions and electrons are in thermodynamic equilib- 
rium, which is not accurate for young supernova remnants. How- 
ever, this is not a major issue for our simulations, as we are inter- 
ested in the large scale, evolved structure of the W50 nebula. A 
comparison of both cooled and uncooled simulations is given in 
Figure |4] 



3.1 Setting the scene for SS433's progenitor 

We make the simple but eff'ective assumption of approximately cos- 
mological abundances of Hydrogen (90%) and Helium (10%) by 
numbeij^ such that = 1.3 is the mean mass per particle in atomic 
mass units. We use a Galactic exponential density profile for the 
ISM, adapted from |Dehnen & Binney 1(1998) ) of the form: 



(4) 



where the constant - 4kpc, = 5.4 kpc is the scale length of 
the stellar disc, and = 40pc is the scale height from the Galaxy 
disc. The pref actor po is determined from a normalisation condition 
that the density at the location of SS 433 is: 



Po -Pism(^« 



ISM V-"'SS433 ' ^SS433 



) = no Jimp 



(5) 



where is the proton mass, and hq is a parameter (for in CGS 
units this defaults to 1 particle per cc, as shown in Figure [TJd). The 
cylindrical coordinate^of SS 433 relative to the Galactic centre are 
approximately R^ = 5.1 kpc and z+ = 204 pc as shown in Figure [T] 

Ordinarily, a gravitational term should be invoked to maintain 
static equilibrium of the background density. In this scenario, the 
gravitational centre of the Milky Way galaxy resides far outside of 
the simulated field of view. To simplify this we allow the back- 
ground temperature profile to have an inverse behaviour to that 
of the background density in order to keep the thermal pressure 
constant, which artificially creates hydrostatic equilibrium in the 
background medium. Preliminary tests were carried out on an un- 
perturbed background to verify the stability of this, and the system 
remained static over a period of 10^ years of simulated time. There- 
fore, over the timescales involved in these simulations, we can 
reliably conclude that subsequent hydrodynamic motion is solely 
driven by the supernova and/or jet interaction. 

The subsequent sections ^ and ^ are concerned with the 
specifics of the supernova and jet models respectively, and the re- 
sults are explored in ^6] iTland ^8] 



^ Courtesy of the University of Chicago Center for Astrophysical Ther- 
monuclear Flashes: http : 1 1 flash.uchicago.edu / website / home / ' 
^ The site of jet-launch in SS 433 is thought to be very much closer to the 
black hole than can be revealed by the resolution of modern telescopes. 
Based on the Schwarzschild radius for a M© black hole, base of the jets in 
SS 433 are likely to be of order 10^^ times closer to the singularity than the 
lobes ofW50. 



^ We note that this corresponds to approximately 70% Hydrogen and 30% 
Helium by mass, which is overabundant in Helium by a few percent with 
respect to the cosmologically accepted primordial values. 
^ This is calculated using <io-GC = 8 kpc as the Sun-Galaxy centre distance, 
<^o-SS433 = 5.5 kpc as the distance from the Sun to SS 433, which has Galac- 
tic coordinates (/ss433, ^SS433) = (39.7°, -2.24°). 
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4 MODELLING THE SUPERNOVA EXPLOSION 

Typically in supernova simulations, a self- similar Sedov-Taylor an- 
alytic solution for a strong blast wave is adopted and used to ini- 
tialise the density, temperature, pressure and velocities of the super- 
nova explosion. This has the advantage that, for a given background 
density and explosion energy, the radius of the advancing shock- 
wave behaves predictably with time, according to the Sedov-Taylor 
( |Sedov| 1959| |Taylor| 1 950) blastwave relation: 

^blastW = a\ (6) 

V Po / 

However, this description is appropriate only during the Sedov 
phase of a SNR blast wave, when the mass swept-up by the super- 
nova shock front is comparable to, or greater than the ejected mass 
from the progenitor. For a sensible progenitor mass ejection and 
an ISM particle {ji - \3) density in the range Hq-QA-I cm~^, the 
Sedov phase occurs when the SNR radius falls in the approximate 
range 10^-10^ pc, which is quite large compared to our resolution 
capability. Thus, we model the supernova in the pre-Sedov free- 
expansion phase, and there are two main advantages to this. 

First, it is difficult to guess the initial blast energy of the su- 
pernova explosion which may have produced W 50. To compensate 
for this unknown we test two values for the initial kinetic energy of 
the gas ejected through the supernova blast. We adopt the notation 
whereby supernova model SNRso has an initial blast kinetic en- 
ergy of E^o - 10^^ ergs and model SNRsi has E51 = 10^^ ergs, and 
together these encompass a reasonable range of supernova cases. 
A purely Sedov-like expansion with thse parameters requires 10^ 
years (for the £^51 case) to 10^ years (for the E50 case) to expand 
to a radius Rwso = ^45 = 45pc (the current size of W50's circular 
shell region, as per Figure [T]d). The free-expansion phase is much 
more rapid initially, but eventually the system hydrodynamically 
establishes a Sedov-like behaviour (explained in This method 
reduces the required time for the SNR radius to reach R^s to some- 
thing much more reasonable, in agreement with observational pre- 
dictions of other old SNRs. The free-expansion SNR reaches a 
maximum radial expansion speed of [1-2]% c and as such a rela- 
tivistic treatment is not critical at this stage. 

Second, the scene of our supernova explosion is set within an 
exponential density background from the Milky Way. As such it 
is beneficial to initialise the supernova explosion as early as possi- 
ble in order to fully investigate the eff'ects of the non-uniform ISM 
upon the SNR as it evolves and sweeps up ambient material. This 
is particularly relevant in order to study the East- West asymmetry 
present in W 50. 

4.1 Implementing the Pre-Sedov Supernova Explosion 

We define /^sedov to be the radius at which a sphere of ISM densit}]^ 
Po contains the mass Mej equal to that ejected from the progenitor, 
and the free-expansion phase is valid for R(t) < /^sedov The ini- 
tial radius Ri of the supernova blast must satisfy two conditions: 
(i) Ri should be as small as possible in order to capture the very 
earliest evolution of the explosion within a non-uniform ISM, (ii) 
Ri must be large enough so that the initial supernova blast ade- 
quately resembles an unpixellated circle on the simulation grid, to 
prevent propagation of unwanted artefacts commonly encountered 



^ The background density can be considered constant across the small 
scales relevant here. 
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when modelling using cartesian grids. We parameterise the initial 
radius as a function of the Sedov radius by 

/3Mei \5 

/?i=/r/?Sedov=/r(^) (7) 

and we find that f = 0.25 satisfies the two conditions above, for 
a mass ejection of Mej = 5Mo and an ambient density of = 1 
particle cm~^. 

The free-fall supernova blastwave is initialised as an over- 
dense discj^of radius Ri about the explosion epicentre (Xc,yc), with 
density 

Pi=/r"'pO. (8) 

All blocks within the radius a^^f = c^efRi of (Xc,yc) are maximally 
refined, so that the grid region in the vicinity of the supernova blast 
is always described by our maximum resolution. Some preliminary 
tests showed c^^f = 1.2 to be sufficient. Using the reasonable ap- 
proximation that the supernova blast energy is initially entirely ki- 
netic, we explore two velocity profiles; the first velocity profile as- 
signs a constant speed Vi to all cells across the supernova initialisa- 
tion zone, with radial unit vector direction: 

More realistically, very young SNR can be modelled as uniformly 
expanding spheres (homologous expansion), whereby (within the 
sphere) the radial velocity of the expanding gas is proportional to 
the radial distance from the sphere centre. Thus, the second veloc- 
ity profile allows the gas speed to increase linearly (a technique also 
used by jJun & Norman| ( [1996| )) from zero at the explosion epicen- 
tre, to a maximum value k^vi at Ri: 

V2{r)r = k/^r (10) 

where = VlO/3 is the necessary constant to maintain the to- 
tal kinetic energy as ^biast- Although preliminary tests indicate that 
these profiles are indistinguishable from one another at large scales 
(see the second velocity profile minimises the effects of in- 
ternal shocks and instabilities to which the supernova explosion is 
susceptible at very early times. As such, the latter velocity profile 
provides a neat method to ensure that the ensuing SNR is more 
closely symmetrical. This is particularly important when consider- 
ing that jets will later be incident upon structures internal to the 
SNR, and will help to minimise stochastic effects which could 
aff'ect the jet propagation in an asymmetric manner, voiding any 
meaningful comparison between the east and west jet evolution. 

We follow the supernova evolution until the blastwave radius 
reaches approximately 44 pc, allowing for some further expansion 
of the SNR towards R^s once the jets are switched on. 

5 MODELLING SS 433'S JETS 

We consider the behaviour of SS433's jet at the current epoch, 
which consists of a jet precessing around a cone. We stress how- 

^ The kinetic energy injected into the SNR disc in the simulation domain 
is a factor ;^;d = (4Ri/3 Sx) of £^biast. due to the ratio of the volume of a 
sphere of radius n with the volume of the disc of thickness Sx used in the 
simulation domain. 
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Jet Models 1-4 












Jet Model 5 snapshot at (l){t) = 
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Figure 2. A graphical summary of the jet models. We use blue and red colours to denote SS 433 's east (mostly blueshifted) and west (mostly redshifted) jets 
respectively, from the point of view of an observer at Earth. This use of colour is purely to differentiate between the east and west jets in the figure, and we 
remind the reader that our simulated jets spend an equal amount of time in approach and recession in the simulation frame, due to the axis of symmetry. The 
green section in (c) indicates the simulation slice to which the precessing jets contribute plasma. In the geometry of (a) and (b), the z axis points out of the 
page, and makes an angle (90 - ^inc) with the line-of-sight towards Earth. 



Table 1. Summary of the jet parameters. 



Symbol 


Description 


Constraint 


Mjet 


Jet mass loss rate 


> 10-6 Mo yr-^ 


Vjet 


Average jet speed 


0.2647 c 


^prec 


Jet cone half-angle 


20.92° 


^prec 


Jet precession period 


162.375 days 


^prec 


Jet angular speed 


3.8695x10-2 radss 


Tjet 


Jet temperature 


10"^ Kelvin 



Numerically, we introduce SS 433 's jets by including relevant 
source terms in the Euler hydrodynamic conservation equations. 
The conservation equations governing SS 433 's jets can be written 
as: 



dp 

dt 



■ V • (pit) 



Pjet 



(11) 



^ In fact, it is much more likely that SS433's jets behaved somewhat diff'er- 
ently in the past, see Section Q 



+ V • (pv(8) v) 



+ V • (p6V) 



■ VP + Pjet Vjet 



■ V-(Pit)-6rad + 6je 



(12) 



(13) 



~W 
dpe 
~W 

where p, v, p, and e are the mass density, velocity, thermal pres- 
sure and total energy density respectively, and 6rad is the radiated 
power (per unit volume) due to radiative cooling. A non-relativistic 
treatment for the jet gas is adopted, i.e. the ISM and jet plasmas 
are assumed to have the same adiabatic index 7 = 5/3. The density, 
temperature, and velocities of the jet injection region are reset to 



(Pje 



Vy ) for each time- step (v^ = everywhere), and the 



ever, that this need not have always been the cas^ and we must 
allow for the possible evolution with time of SS433's jet preces- 
sion since the formation of the compact object. We discuss various 
interesting jet geometries in turn, in the subsections that follow. 
The parameters which govern jet behaviour in our models are sum- 
marised in Table [T] along with their default values in convenient 
units (note that SI units are assumed wherever these parameters ap- 
pear in equation form). 

Note that Mjet and Pjet are poorly constrained. Since the mass 
loss in the jets is not well known, we investigate three orders of 
magnitude in Mjet ranging from a relatively weak jet IQ-^ M© yr.-\ 
an intermediate jet of IQ-^ M© yr-\ and a strong jet of IQ-"^ 
Mq yr.-^ The jet temperature Pjet is also not very well known, but 
the prominent H^^ emission in the jets suggest the temperature is at 
le ast 10"^ Kelvin. An estimate of SS 433's jet temperature is given 
in |Fabrika| ( [2004t as 2 x 10^ K. 



density pjet is a model-dependant function of Mjet- 

The structure created by a jet in its surrounding environment 
is primarily determined by the jet velocity components Vx, Vy and v^ 



which are also model dependant, but their most generalised forms 
follow the kinematic modeOof SS 433 such that: 



Vx(^jet(0, 0(O) = Vjet COS (^jet(0 + Q) X 

Vy[e-^a{t), (p(t)) = Vjet sin (^jet(0 + Ce) cos 0(0 y 
v,(^jet(0, 0(O) = Vjet sin (6'jet(0 + Q) sln 0(0 I 

0(0 = <^prec (t - to) + 00 

( for the Western jet 
7T for the Eastern jet 



(14) 

(15) 

(16) 
(17) 

(18) 



where 0(0 is the precession phase angle at time t, and is the time 
when the jets switched on. Throughout this paper we use the defi- 
nition whereby precessional phase zero occurs when the red-most 
(western) jet is maximall}^] red- shifted (i.e. the blue-most jet is 

Note that here the kinematic model is adapted slightly because SS 433's 
mean-jet-axis is the x-axis of the simulation domain. This explains the ab- 
sence of the inclination angle dependence cos OiiYc 

and sin Oinc ill this repre- 
sentation. 

^ ^ In the simulations, the maximum redshifts are equal for the east and west 
jets, due to the axis of symmetry used. 
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maximally blue- shifted). SS 433 precesses in a clockwise fashiorp] 
when looking along the direction of the western jet from the view- 
point of SS433. Consequently, if the phase angle starts from zero 
diit - to then the offset 0o = -njl ensures that 0(0 = nnjl when 
the jets are maximally Doppler shifted, and 0(0 = nn when the jets 
display only the transverse Doppler shift (in the xy-plane), where 
{n G Z*} is the revolution number. 

To allow for the possibility that the precessional motion in 
SS 433 has changed during its lifetime as an X-ray binary, either 
smoothly or as a discontinuous process (c.f. episodic jet outbursts 
in Cyg X-3 mentioned in g2.3| ), we introduce into our models some 
simple time-dependence in the precession cone angle: 



cjet 



(0 = ^0 + (^ - ^o) ^ 



(19) 



which can be made constant by setting ^=0. 

The jet density depends on the jet mass loss rate Mjet and some 
model-dependant geometry, and the internal pressure is calculated 
from the density: 



Pjet ~ Q^geom -^jet 

_ Pjet ^jet 
iet — 



(20) 
(21) 



where ofgeom has units of time I volume. The total power per unit vol- 
ume for the jets has kinetic and thermal components: 



(r-i) 



Pjet( 



B ^jet 



Idmpiy- 1) 



■)• 



(22) 



We explore 5 different jet initialisation models, and later dis- 
cuss the relative successes of each. Models 1 to 4 are based upon the 
jet motion described in equations ([14]) through ([T9| and depicted in 
Figure [2^, whereas model 5 attempts to mimic SS433's discrete 
ballistic jet ejecta, as depicted in Figure [2|3. 



5.1 Jet Model 1 - Static Cylindrical Jet 

Our most basic jet model consists of a static (non-precessing) cylin- 
drical jet along the x-axis. This is a special case of the general equa- 
tions ([14]) through ([19J, created by setting ^0 = 0, 6* = 0, and 0(0 = 0, 
such that the velocity components are reduced to = vjet cos Cq x 
and Vy = = 0. 

The density of the jet is simply a product of three factors: 
the mass loss rate Mjet in the jets, the (resolution-dependant) pixel 
crossing time - (5x/vjet, and the model-dependant jet volume 
(;ijx Wjy ^jz Sx^), where {n-^^, «jy, ;ijz) are the number of pixels consti- 
tuting the jet initialisation region in the respective (x,j,z) direc- 
tions. This issue can be simplified by choosing the initial jet size 
(as in Figure|2^) to be 1 cubic pixel {n-^^ -n-^y-n-^^^-V). This makes 
the most physical sense, since the pixel size is very much larger 
than the jet structure revealed by VLBI imaging, which is in turn 
very much larger than a few Schwarzschild radii of the black hole 
in SS433. Thus, the density in this model is given according to 
equation ([20| with: 



1 



2 Vjet 



(23) 



The right hand grip rule is useful here: gripping the x-axis with the right 
hand, with the thumb in the +x direction, the curled fingers indicate the 
precession direction for the west jet. 



where the above summed over both jets will give the required mass 
loss rate per unit length (one pixel) in x. 

Simulations based upon this jet model were performed for 
each of the three Mjet values, both with and without the presence of 
an evolved SNR. The simulation runs for this model are also used 
to quantify the east- west jet asymmetry in W50 by constraining 
the density normalisation uq and Galaxy scale-height Zd in equa- 
tion ([4]) appropriate at the location of SS 433 in our Galaxy, and 
these results are discussed in g7.1| 

5.2 Jet Model 2 - Conical Phase-shiftng Jet: ( 6-^^ = ) 

This model invokes a time-averaged representation of SS 433 's pre- 
cessing conical jet, termed phase-shifting. As with model 1, the 
east and west jets are each just 1 cubic pixel in size. This model 
uses present-day observed parameters describing SS433's preces- 
sion at the current epoch, and assumes the jets have not changed 
throughout their history. 

The computational time-step is adjusted so that precession 
phase-space is sampled times per precession period: 



A^st 



A0. 



*^prec 

In 



step 



(24) 



(25) 



If we define a phase-tolerance as 0toi=A0step/2, then the precessing 
jet contributej^to the simulation domain when the phase 0(0 falls 
within +0toi of (for the upper-half of the plane) or n (for the lower- 
half). When this happens, the phase 0(0 is rounded to the nearest 
of either of (0, n), and this is equivalent to introducing an error of 
+(50/;is) percent to the precession period (hence the term phase- 
shifting). 

The velocities are initialised by setting ^0 = ^prec, ^ = and 
0(0 = <^prec^ in equations (T4| through ([19]), where co^^qc - ^^/tproc- 
Assuming that the amount of mass outflowing in the jets in one 
precession period can be written as Mjet ^prec, the density is set ac- 
cording to equation ([20|) with: 



(26) 



2 n. 6x^ 



We choose n^ = 20 as a compromise value to prevent impracti- 
cal inflation of the computational time, whilst limiting the error in 
precessional period to just +2.5%. 



5.3 Jet Model 3 - Conical Phase-shiftng Jet: ( djet > ) 

This model is similar to Jet Model 2 in that it follows the same 
precession phase -shifting model, but with the exception that the 
precession cone angle is allowed to increase linearly with time. 
The basis for this model comes from the discrepancy between the 
opening angle of the jet precession cone ^prec currently observed in 
SS 433, and the angle subtended by the lobes of W 50 at the central 
source SS 433. The green lines in Figure[T]3 represent the current jet 
precession angle projected out to the extent of W50. Had SS 433's 
jets penetrated the lobes of W50 at this angle, we would expect 
quite a different morphology to result. To investigate the effects of 
a smoothly evolving jet using equation |T9] l, the initial half-cone 
angle is set to ^0 = and the derivative becomes 



The jet pixels are set according to the jet values (pjet, Tjet . ^jet) rather than 
the background ISM values. 
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(27) 



such that the cone half-angle increases linearly to the currently ob- 
served value, where /jets is the time taken for the jets to reach the 
extent of the lobes in W 50. The jet density is as described in model 
2. 



5.4 Jet Model 4 - Conical Phase-shiftng Episodic Jet 

In this model we investigate the possibility that SS 433 's jet activ- 
ity is discontinuous. This is inspired by the microquasar Cyg X-3, 
which exhibits episodic jet outbursts, each with quite different jet 
speeds and precession cone angles. To emulate this behaviour, we 
use the same phase-shifting model and jet density as described in 
jet model 2, but the outbursts are created by using ^ = in equation 
( p^ and setting according to: 





for 


ti <t <t2 






for 


t3 <t <t4 


(28) 


03 


for 


ts<t< te. 





5.5 Jet Model 5 - Fireball Jet Model 

The fireball model was inspired by the spectacular animatioip] 
showing the discrete fireball-like nature of SS433's jets from the 
work of |Mioduszewski et al.|(|2003| ). The optical spectroscopy ob- 
servations of Blundell et al. estimate the expansion speed Uq 
of the jet fireballs in the innermost regions of SS 433 as approxi- 
mately We - 1% of the jet speed. Using this estimate of the fireball 
expansion rate, and by assuming that the initial fireball size upon 
ejection is negligible compared to our resolution, we approximate 
the initial jet collimation angle SLSif/ ^ We/vjet - 0.01 radians. Thus, 
it is possible to calculate the distance from the origin at which the 
fireballs will have expanded to a size comparable to 1 simulation 
pixel, such that we can reasonably resolve the discrete ejecta seen 
in VLBI observations on the simulation grid. For approximately 
spherical fireballs of radius i?fb then (2Rfb = rfb = 6x) gives the 
appropriate fireball displacement rfb along its trajectory r. Hence 
we can write the coordinates of the fireball pixel in each quadrant 
as: 

6x cos ^nrec 

Xfk = xo ± — = xo ± 93A6X (29) 



yfb 



yo ± 



Sx sin 



- yo ± 35.7 Sx 



(30) 



and there are 4 such fireball pixels; one in each quadrant of the 
x_y-plane as indicated in Figure |2]3. 

In three dimensions the jet sweeps out part of a corkscrew 
structure in space over one precession period, however our domain 
is eff'ectively a slice one pixel (6x) thick in the z direction. It is use- 
ful to define the quantity /fb which is the angle subtended by one 
fireball pixel about the x-axis, as a fraction of one whole preces- 
sional period (In radians): 



fa. 



Sx 
Inrfk 



In sin Or,, 



(31) 



where 6x is the pixel size and eff'ective width of the simulation slice 
in z. We now introduce a user-defined constant /Cfb (equal to 1 by 



See: http : / /www. nrao.edu/pr/ 2004/ ss433/ ss433.movie.gif 



default), which allows the user to establish a compromise between 
computation time and resolution of the fireballs. The precessing 
jet only contributes plasma when it sweeps through the simula- 
tion domain slice. In phase-space this spans half of the phase- width 
0w = 27TKfb /fb each side of the x_y-plane, or (0 + 7r/Cfb/fb) for the fire- 
balls in the northern (upper-half) plane, and (tt ± 7r/Cfb/fb) for the 
southern (lower-half) plane. This phase range is specified by the 
criterion: 



cos^ 0(0 > cos^ 



(32) 



The computational time- step A^step must be reduced for this 
model, to ensure that the jet phase 0(0 evaluated at and + A^step 
does not cause the jet to skip the x_y-plane altogether. The preces- 
sion crossing time for the jet to sweep through the simulation slice 
is equal to A^cross = ^fb /fb ^prec- To ensure adequate sampling of the 
precession period, the computational time-step in this model needs 
to be equal to or smaller than the precession crossing time: 



A/^step ~ A/^ej-Qss 



^fb ^pre 

In sin 



; 0.72 days 



(33) 



which (with /Cfb = 1) is typically a few orders of magnitude smaller 
than other dynamical timescales on the simulation grid. The jet 
plasma density is given according to equation ( [20| with: 

^prec 



26x^ 



(34) 



The jet plasma is initialised for the fireball pixels according to: 

(Pjet , Tjet , Vjet ) for COS^ 0(0 > COS^ 0^ 



(piSM , T^iSM , O) Otherwise. 



(35) 



If the criterion of equation ( [32] ) is met (thus signifying a jet 
intersection with the xj-plane), the precessional phase must be very 
close to either of (0, tt) and we round 0(0 to the nearest of these. 
This ensures that tiny residual values are set to zero when the 
fireball velocities are initialised: 



y(0(O) 
z(0(O) j 



COS(^prec + Cq) 

sin(^prec + Ce) cos 0(0 
sin(^prec + C^) sin 0(0 , 



(36) 



where Cg has the usual meaning of zero for the west jet (x > 0) and 
TT for the east jet (x < 0). 



5.6 Computation and Data Analysis 

A total of 32 preliminary simulations were conducted for this pa- 
per, all of which were performed using the Oxford Astrophysics 
Glamdring cluster. The characteristics of each simulation are out- 
lined in Tables|2]and|3] and for convenience we will be refer to each 
simulation by its shortname given therein. 

All of the scripts used to analyse the data from our simula- 
tions were coded using the Perl Data Languag^ The simulations 
employ the Hierarchical Data Format (HDF5) option for data out- 
put in FLASH, which is appropriate for use with AMR. Using code 
written in PDL, each HDF5 file is converted to a uniform grid for 
easier analysis, and stored as (archived) Flexible Image Transport 

The Perl Data Language (PDL) has been developed by K. Glazebrook, 
J. Brinchmann, J. Cerney, C. DeForest, D. Hunt, T. Jenness, T. Luka, R. 
Schwebel, and C. Soeller and can be obtained from http://pdl.perl.org. 
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Table 2. Summary of our preliminary simulations testing the SNR parameters. The first six runs test the background profiles and radiative cooling. The 
last four simulation runs are variations of SNR6, where meaningful astrophysical parameters are tested. 



Shortname 


Description 


Parameter Details 




SNRl 


Supernova only 


SNR defaults and: p-profile = 0, v-prof = 0, Cooling: 


=off 


SNR2 


Supernova only 


SNR defaults and: p-profile = 1, v-prof = 0, Cooling: 


=off 


SNR3 


Supernova only 


SNR defaults and: p-profile = 0, v-prof = 1, Cooling: 


=off 


SNR4 


Supernova only 


SNR defaults and: p-profile = 1, v-prof = 1, Cooling: 


=off 


SNR5 


Supernova only 


SNR defaults and: p-profile = 0, v-prof = 1, Cooling: 


=on 


SNR6 


Supernova only 


SNR defaults and: p-profile = 1, v-prof = 1, Cooling: 


=on 


SNR6b 


Supernova only 


As SNR6 but: £:biast = 10^^ ergs 




SNR6c 


Supernova only 


As SNR6 but: = 0.1 cm'^ 




SNR6d 


Supernova only 


AsSNR6but: Mgj = lOM© 




SNR6e 


Supernova only 


As SNR6 but: Tbg = 10^ Kelvin 





SNR default parameters: Mgj = 5Mo, = 1 cm ^, £^biast = 10^^ ergs, Tbg = 10^ Kelvin, Zd = 40 pc 



System (FITS) files with appropriate header information relevant 
to each simulation. 



6 SNR: RESULTS & DISCUSSION 

A series of 10 preliminary supernova tests were performed to in- 
vestigate the effects of the key physical parameters governing the 
evolution of the SNR (see Table |2]). The supernova blast energy 
£^biast» the background density hq at SS 433, and the mass ejected by 
the progenitor Mej were among the parameters tested, as well as the 
effects of radiative cooling and changing the background density 
profile from uniform to exponential. The velocity profiles described 
in section g4.1| are also tested, but they do not have a meaningful 
interpretation. 

A sample of snapshot images from two of these SNR simula- 
tions, SNR 6c (left column) and SNR 6 (right column) are shown at 
various time intervals in Figure [3] The figure illustrates the differ- 
ence in the SNR evolution when the background density at SS 433 
is varied from no = 0.1 cm~^ (left column) to no = lcm~^ (right 
column). The difference in the initial radii of the supernova blasts 
from the first image in each column is a consequence of the first 
density difference by a factor of 10 in equation ([T]). The most ob- 
vious difference is in the time taken for the SNR to expand to /?45, 
whereby increasing the density by an order of magnitude causes 
the duration of expansion to increase by more than a factor of 2. As 
the supernova Shockwave expands into the ISM, the effects of the 
exponential density gradient become apparent in a number of ways. 
First, the density of the west side of the SNR (nearest the Galaxy 
plane) is noticeably higher than the density at the east side in the 
evolved shell. The ratio of the densities Rp = Pw/Pe on opposite 
sides of the SNR shell reaches a maximum of Rp ^ 10 along the 
direction normal to the Galactic plane (vector GD in Figure [!} in 
the snapshots at later times, in Figure|3] As expected, the east-west 
density ratio has approximately the same value for both SNR 6 and 
SNR 6c. Second, the position of the centre {xo,yo) of the SNR quite 
noticeably shifts as a function of time, with respect to the initial 
blast epicentre (Xc, jc), moving away from the Galaxy plane, down- 
stream in the density gradient. 



6.1 Monitoring the SNR shock front 

To quantify the differences in morphology of the resultant SNR 
produced by each simulation in Table |2] a shock-front detection 
algorithm was applied to each dataset. The coordinates {x\,y\) of 
the points lying on the shock-front were fit to a general ellipse of 
the parametric form: 

Xi(^i) - xo + a cos(0) cos(^i) - b sin(0) sin(^i) (37) 

y.{e-;) ^yo + a sin(0) cos(^i) + b cos(0) sin(^i) (38) 

where {xo,yo) is the centre, is the angle each point (xi,_yi) sub- 
tends between the x axis and the centre in the anti-clockwise sense, 
is the angle between the x axis and the semi-major axis, and a 
and b the semi-major and semi-minor axes respectively. The eccen- 
tricity as well as the displacements of the SNR shell centre (xq, jo) 
from the explosion epicentre (Xc, jc) are then easily calculated as: 

e = Vl - 0^1 (39) 

Ax = Xc - xo (40) 

Aj = jc-Jo- (41) 

The physical variables (density, pressure, velocities etc) are 
sampled at the detected shock-front, and various quantities are plot- 
ted as a function of time or SNR radius in Figure |4] The shift in 
centre position of the SNR are plotted in Figure |4jg) and (h), and 
it is clear from Figure |4|i) that although the SNR becomes slightly 
elliptical in the presence of an exponential density background, the 
effect is a minor one {e < 0.2) compared to the shift of the SNR 
centres. In Figure [ija) the average radius R - Va^ + b^ of the SNR 
is plotted as a function of its age, and the dotted lines show the best 
fit to a f-l^ Sedov-like function of the form k + a[^^ where 
k is some offset in parsec. The shock generated through the super- 
nova explosion is strong as evident from the compressior[^ratio of 
Pilpo ^4 across the shock (Fig.|4]D) as well as the ram pressure be- 
ing higher than the thermal pressure (Fig.|4p, and the shock speed 

Note that the shock-front densities plotted in Figure [i] are actually av- 
eraged values from an annulus 1 pixel thick about the detected location of 
the shock-front, and that densities on the east and west side of the SNR dif- 
fer by an order of magnitude when the background follows the exponential 
density profile. 
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Figure 3. The evolution of 10^^ erg supernova explosions in the presence of the Galaxy density gradient. Supernova simulations SNR6 and SNR6c 
are shown at 4 time intervals during the SNR evolution, plus the final snapshot when the SNR reaches ~45pc in radius. The x symbol indicates 
the explosion epicentre on each snapshot, and the displacements Ax and A_y of the SNR from the epicentre (xc,yc) are also given. RTI indicates the 
occurrence of Rayleigh-Taylor instabilities, where lower-density gas from inner radii are flowing into regions of higher density gas at larger radii 
within the SNR. The density for the left-hand column (SNR 6c) is 10 times lower than the colour bar indicates. The animations associated with these 
simulations are available at: http : / /www-astro. physics.ox.ac.uk/~ptg /RES EARCH/ re search.html 
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Figure 4. A shock-front detection algorithm is used to monitor the size and shape of the SNR, and to sample the density, pressure and velocities 
of the SNR shell as a function of time. The figure key summarises the SNR simulation details given Table [2] where "E" precedes the log of 
the blast energy in ergs, "v" precedes the the velocity profile used (see ^4.1} , "rho" precedes the background density profile used (O=uniform, 
1= exponential), and "cool" indicates whether cooling is on (=1) or off (=0). For simulations SNR6b-SNR6e, parameters are only specified 
where they differ from those used in simulation SNR 6. Note that SNR 6c in a lower density environment reaches 45 pc much faster than the 
other SNRs, and is truncated when the x-axis is time. SNR 6b has an order of magnitude lower blast energy than the others, and its shock front 
dissociates and loses identity before it reaches 30 pc, such that the shock front cannot be well detected and these data are also truncated. 
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remains hypersonic (Mach»5) for the duration of the simulations 
as shown in Fig.|4]: and|4ji. In Figure]^, the mass contained within 
an SNR sphere is roughly estimated by averaging the density of the 
SNR disc of radius R in the simulation domain, and scaling up by 
the volume of a sphere of the same radius. An estimate of the mass 
expected to be swept up by a spherical SNR in a uniform back- 
ground of density uq is also plotted onj^ for comparison, with the 
functional form 

/(Mej,^o,r) = Mej+ "^^P (r'-R^l for r > R, (42) 

such that the initial mass contained within the young SNR of ra- 
dius Ri is equal to the mass ejected Mej by the progenitor. This is a 
crude model, as it makes no attempt to account for the exponential 
background density or the ISM material contained within R^, but 
it proves to be consistent with the scaled-up mass from the SNR 
to within a factor of a few. The final plot (Figure |4]<:) tracks the 
total energy of the system, which is again scaled up by a factor 
- 4/?i/(3 6x), representing the ratio of the volume of a sphere 
of radius R^ to the volume of the initial disc of the same radius, 
used to initialise the supernova explosion in the simulation domain 
(see ^4.1 1 for details about this). The projected energy in a sphere 
of radius /?snr(0 as shown in Figure |4]<: is then: 



(43) 



representing the discrete summation of the internal and kinetic en- 
ergies per unit volume, for concentric annuli of volume elements 
dVi from the centre of the SNR to /?snr(0- The energy radiated via 
cooling is tracked as a function of time for each dataset, and radia- 
tive loss of energy from the system via cooling is negligiblj^for 
all but one of the SNR simulations. SNR 6b, with an order of mag- 
nitude lower blast energy than the others, was found to slow down 
much faster than the other SNRs. This is expected for a remnant 
passing through the Sedov phase into the radiative phase in which 
cooling does become important, and the addition of cumulative ra- 
diative losses is included in the total energy budget of the system, 
as indicated by the dashed blue- white line in Figure [4JC. SNR 6b 
eventually began to disperse and lose its circular SNR identity, at 
which point the shock-front detection algorithm failed to detect a 
well-defined shock and for this reason the blue lines in Figure|4]are 
truncated at 28 pc. The other 9 SNR simulations can be thought of 
as adiabatic to a good approximation, evidenced by the constancy 
of the total energy £tot shown by the solid lines in Figure[4]<:. As one 
might expect, the shock deceleration is less noticeable in the lower 
density background of SNR 6c, and the transfer of kinetic energy 
into thermal energy is more gradual, as Figure [4]<: shows. 



6.2 The displacement of SS 433 from the SNR centre 

The shockfront detection algorithm applied to each of the SNR sim- 
ulations allowed us the fit an ellipse to the shockfront and track the 
movement of its centre away from the blast epicentre, as a function 
of time (Figure |4^ and h). 

The movement of a SNR's centre in the presence of a den- 
sity gradient was mentioned by [Chevalier & Gardner | { [1974| ). The 
cause of the shifting is not buoyancy, but rather the relative ease 
with which the eastern side of the SNR can propagate in the lower 



Note that the microphysics of emission-line cooling is not yet fully im- 
plemented in this version of the code. 



density ISM, compared to the much higher inertial resistance met 
by the western side of the SNR in the higher density towards the 
Galactic plane. In other words, the shifting is an apparent eff'ect due 
to the independent propagation of the east and west shock fronts 
into diff'erent media, and the SNR shell doesn't physically move as 
a rigid body. 

According to SNR 6 of Figure by the time the SNR ex- 
pands to 40-45 pc in radius we can expect the centre position to 
have shifted by {^S^,^Sy) ^ (-7.1,-2.5) pc or (-4.4,-1.6) ar- 
cmin. Note that these shifts are taken from Figure |3] and therefore 
are only appropriate to the SNR evolution in an environment where 
the galactic scale height is Zd = 40 pc, and the density normali- 
sation is fiQ - 1.0 cm~^. The displacement values for the case of 
SNR 6c diff'er only slightly from those of SNR 6, but the displace- 
ment is likely to be larger in the case of a steeper background den- 
sity gradient (e.g. Zd = 30 pc). Regardless of the exact value of 
Zd, the apparent shift in SNR-centre will always be in roughly the 
same direction with respect to the density gradient, i.e. away from 
the Galactic disc with -ve A5x and A^Sy. 

We apply a rotation matrix through the clockwise angle 6a 
from Figure [TJ) to determine the equivalent SNR shifts along the 
ordinates of right ascension and declination: 



A5, 



cos 9a sin 9a 


i -AS, \ 


' +4.6 ^ 


- sin 9a cos 9a 


[ AS, J 


. -1-0 . 



arcmin(44) 



where the -A^x occurs because the simulation x-axis and the right 
ascension coordinate are almost antiparallel. 

Without understanding that the SNR shell hydrodynamically 
shifts position in the density gradient of the Galaxy, one might mis- 
takenly calculate that SS 433 is moving away from the SNR cen- 
tre. The movement of the SNR must be considered when calcu- 
lating the age of the black hole by a simple means of off'-centre 
displacement from the SNR centre, divided by SS433's apparent 
proper motio n. We refer the reader to the detailed study of W 50 

( |2007|^' 



described by 



Lockman et al. 



where the off'-centre dis- 



placement of determined from X-ray observations ( [Watson et al.| 
|1983| ) is used in conjunction measurements of SS 433 's proper mo- 
tion (A. Mioduszewski & M. Rupen, private communication), to es- 
timate the age of W 50 of order 10^ years. 

The correct way to estimate the age of W50 by this 
method must include both the proper motion data from SS 433 and 
a contribution from the apparent SNR displacement, such as: 



^S a - i/J^a -^a) %50 
A^'^ = (lis - Vs)tw50 



(45) 
(46) 



where jia and /us describe SS 433 's proper motion, Va and vs are the 
time-averaged SNR displacement speeds from the simulations ac- 
cording to Va = ASa/tsNR and Vs = ASs/tsNR, and ^snr is the appro- 
priate simulation time over which the displacement occured. Equa- 
tions ( [45| l and \46\ should give the same answer when solved for 
%50, but this method requires an accurate measurement of AS a and 
A^^, which have not yet been obtained from observations. 



The displacement values are not used here because the article quotes the 
linear or angular displacement as 4 pc or 5 arcmin. However, these length 
and angle values are discrepant by a factor of 2 for objects at a distance of 
5.5 kpc. 
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Table 3. Summary of our preliminary simulations testing the jet parameters. Simulation Jetl is used as the control group to which Jet2-Jet8 are compared. 
Jet9 investigates the effects of an episodic (rather than continuous) jet, and Jet 10 demonstrates the fireball model described in the text. The Galaxy scale 
height parameter Zd is tested in another series of simulations, as detailed in Figure[6] 



Shortname 


Description 


Parameter Details 


Jetl 


Jets only (Model 1) 


Jet defaults and: ^0 = 0°, 6* = 0° yr'^ 


Jet2 


Jets only (Model 1) 


As Jetl but: Mjet = 10"^ Moyr"! 


Jet3 


Jets only (Model 1) 


As Jetl but: Mjet = lO'^Moyr-^ 


Jet4 


Jets only (Model 1) 


As Jetl but: no = 0.1 cm~^ 


Jet5 


Jets only (Model 2) 


As Jetl but: = 10° 


Jet6 


Jets only (Model 2) 


As Jetl but: = 20° 


Jet7 


Jets only (Model 2) 


As Jetl but: = 40° 


Jets 


Jets only (Model 3) 


As Jetl but: 6* = 20°/1000};r 






6*1 =40°for0<r< 1000 yrs 


Jet9 


Jets only (Model 4) Jet defaults but = 0.1 cm~^ and: 62 = 0° for 1000 < t < 2600 yrs 






O3 = Oprcc for 2600 < t < 2700 yrs 


JetlO 


Jets only (Model 5) 


Jet defaults but no = OA cm~^ and Kfk = 10 




Jet default parameters: Mjet 


= 10-4 Mo yr-\ no = 1 cm-^, Zd = 40 pc 



7 JETS: RESULTS & DISCUSSION 

A series of 10 simulations were performed to investigate the effects 
of the jet kinematics, such as the jet precession cone angle and jet 
power (influenced through Mjet), and the details of the simulations 
are given in Table [3] Each jet simulation was allowed to continue 
running until it became obvious that a particular jet model was un- 
likely to produce the required morphology observed in W 50. An 
image showing the density map for each simulation in its final stage 
is shown in Figure[5] and a contour map of the |Dubner et al.|(T998| ) 
observation of W 50 has been overlaid upon each density map for 
comparison. The magnitude of the jet velocity is kept constant at 
Vjet = 0.2647c for each jet simulation. 

Figure [5^ shows Jet 1 from Table [3] consists of a simple static 
cylindrical jet along the x axis, with a mass loss rate of Mjet = 
lO-'^Moyr"^ in an exponential Galactic density gradient of scale- 
height Zd = 40 pc and normalised to /iq = 1 cm"^ at SS433's lo- 
cation on the grid. Although the east- west extent of the lobes from 
Jetl approximately matches that of W50, the eastern jet lobe has 
reached a displacement of ~125 pc from SS 433 which is in excess 
of W50's eastern lobe extent, whilst the simulated western jet lobe 
is short of W 50' s western lobe, as clearly indicated by the contour 
line. This could indicate that the galaxy scale-height used was too 
small, and this is investigated in ^7.1| As the jets plough through 
the ISM, they carve out a cavity of lower density which enables 
the easier propagation of the proceeding jet material. Referring to 
the colour scale used in Figure |5^, a contrast of an order of mag- 
nitude or more in density is achieved in some places, between the 
lower density ISM directly surrounding the eastern jet (blue) and 
the denser gas of the cavity walls (red) behind the jet shock front. 
We can define an average lobe speed as the ratio of the average ex- 
tent of the lobes to the time taken for the jet shock front to reach 
such an extent: 



where viobes = 0.148c for Jet 1, based on 200 pc as the east- west 
lobe extent and 2200 yrs of simulation time. Note however, that 
since the jet is "on" for the duration for this simulation, that the ki- 



netic energy of the jets is continually replenished, and the preceding 
jet activity has already evacuated a cavity in the ISM. This results in 
the speed of the jet itself remaining very close to Vjet until very near 
to the bow shock at the ends of the jet lobes, and this is described 
further in |8.7| Upon reaching the ends of the lobes, the jet ejecta 
decelerates rapidly as momentum and kinetic energy are dissipated 
into the denser ambient ISM, to drive the shock forward. As the jet 
shock-front propagates, it imparts momentum into the transverse 
(y) direction, such that the wake of the jet has a well-defined width 
at a given distance from the jet launch point. The result is that even 
a cylindrical jet (zero cone angle) has a substantial width as evident 
in Figure |5^, where clearly the western jet lobe is as wide as (and 
in some parts, too wide for) W50's western lobe. The width of the 
jet lobe as a function of distance from SS 433 is a useful diagnostic 
in constraining the nature of the jet which sculpted W 50's lobes, as 
shown in Figure [To| 

Jet 2 of Figure |5]3 is identical to Jetl but for the mass-loss 
rate which is an order of magnitude lower at Mjet = 10~^ M© yr~^ 
The level of jet collimation displayed by Jetl is not true for Jet 2, 
and although the jet mass loss rate is still relatively high it is clear 
that, if allowed to expand to double the size shown in Figure |5]3, 
the peanut-like lobes created by Jet 2 would be much too wide to 
produce the lobes of W 50. In this case the jets have significantly 
decelerated before moving far from the jet launch point, and the 
jet loses its well defined shape, bending and buckling and thereby 
transferring relatively more momentum in the transverse direction 
than does Jetl. The averaged lobe speed for Jet 2 is viobes = 0.058c. 
As one might expect, the deceleration and lack of collimation are 
both even more obvious in Jet 3 of Figure [5j:, where the jet mass- 
loss rate has been reduced by a further order of magnitude to Mjet = 
10~^Moyr~\ resulting in an elliptical bubble with average speed 

Viobes = 0.017c. 

Jet 4 of Figure |5]i has the same jet mass loss rate as Jetl, but 
in this case the ambient density normalisation is lower by an order 
of magnitude: no = OA cm~^. The relative ease with which this jet 
propagates through the ISM is apparent both through the level of 
collimation of the jets and also the jet lobe extents exceed those 
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Figure 5. Final snapshots are shown of a series of jet simulations demonstrating the effects of relevant jet parameters upon the morphology of the jet cocoon. 
Refer to Table [3] for a summary of the differences between these simulations, and ^for a detailed commentary on the interpretation of each one. The 
animations associated with these simulations are available dLi: \http : 1 1 www-astro. physics. ox. ac.uk I ~ptg / RES EARCH/ research.html 
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of W 50 on each side. As one might expect, the jets require only a 
relatively short time (less than 1650 years) to propagate to the full 
extent of W50's lobes (viobes = 0.219c), as the background den- 
sity is lowered and we tend towards the in vacuo limit described 
in |2.2| The transverse momentum component is reduced relative 
to the axial component, and jet travel time is short enough that the 
wake of the simulated jet lobes are thinner than the lobes of W 50. 
This leads to an interesting relationship between the thickness of 
the jet lobes and the ratio Mj^Jno. The effect upon the jet lobes 
morphology of reducing no by a factor of 10 is equivalent to in- 
creasing Mjet by the same factor. For example, repeating Jet 4 with 
Mjet reduced by a factor of 10, reproduces the same jet morphol- 
ogy as Jet 1, and running Jet 1 with an order of magnitude decrease 
in the background density, reproduces the same jet morphology of 
Jet 4. This relationship is described further in ^7.1| 

Jet 5 of Figure [5^ utilises Jet Model 2, and has similar param- 
eters to Jet 1 with the addition of time-averaged precession of the 
jet with half-cone angle 10°, approximately half the jet cone angle 
^prec of SS 433 current precession state. Even with this relatively 
small half-cone angle it is easy to see that the evolved lobes would 
be much wider than W50's lobes, but despite the peanut shaped 
morphology this jet still produces two distinguishable lobes. Inter- 
estingly Jet 6 of Figure|5f features a rounded central region between 
two stubby lobes, when the jet half cone angle is set to ^prec- When 
the precession half-cone angle is increased to 40° is with Jet 7 of 
Figure [5^, the lobe structure is lost and the jet cavity more closely 
resembles a circular blown bubble. Each of the simulations with 
non-zero precession cone angles display evidence of hydrodynamic 
refocusing of the jets towards the mean jet axis (for more details see 

To investigate the possibility that SS433's precession devel- 
oped gradually over time. Jet 8 of Figure [5}i features a linearly 
growing jet with half-cone angle starting from zero, and reaching 
^prec after 1000 years. It was thought that this jet in early stages 
(with a small precession cone angle) could produce the desired 
elongated lobes of W50, and then at later times when the preces- 
sion cone angle is larger, the jet would assist in inflating the circular 
SNR-like region. Interestingly, the result is quite different because 
the jet interacts with itself. As before, the jet ploughs through the 
ISM and carves out a cavity of low density, providing much lower 
resistance for the jet that follows at later times. However, the dense 
shock-front of ISM gas swept-up by the early jet, also acts to con- 
fine the later jet. Rather than penetrating the cavity walls, the jets 
interact with the shell at a grazing angle, and are "guided" back to- 
wards a path of lower resistance. More details of this effect follow 
in^ 

7.1 Jet lobe propagation: constraining SS 433's environment 

The east- west asymmetry in W50's jet lobes, if due to the Galaxy 
density gradient, provides us with an opportunity to do two things: 
(i) to investigate hydrodynamically the large-scale effects of jet 
propagation in a non-uniform background density environment for 
comparison with observations, and (ii) to constrain the parameters 
of the Galactic density profile appropriate to the location of SS 433. 

The purpose of this section is to constrain the parameters no 
and Zd by comparing the jet lobe extents from the simulations to 
those observed in W 50. These two parameters have the following 
effects: 

(i) Density normalisation no - this defines the relative ease with 
which the jets propagate through the ISM. Increasing no causes 



120 



100 

o 

^ 80 

B 
>< 
w 

CD 



20 







Figure 6. We monitor the east and west lobe lengths as a function 
of time, whilst varying the Galaxy scale-height Zd and ambient den- 
sity no. The distinctive morphology of W50 can thus be used to con- 
strain the physical parameters most appropriate to the SS433-W50 sys- 
tem. The jet shock-front detection technique can be seen in action at: 
\http : 1 1 www-astro. physic s. ox.ac.uk I ~ptg IRES E ARCH I re search.html\ 

greater deceleration of the jets, because both the east and west 
jets experience greater resistance from the ISM through which they 
propagate. The momentum lost from the jets during deceleration 
is transferred to the ISM, with a large portion of the momentum 
transfer being in the direction perpendicular to the mean jet axis. 
As a result, the jets become lesser collimated, and in extreme cases 
the jets create a peanut-like cocoon shape, rather than well defined 
lobes. 

(ii) Galaxy disc density scale-height Zd - this scale height de- 
termines how fast the density decreases with distance from the 
Galaxy plane, such that moving further away from the Galaxy plane 
by Zd parsec causes the density to change by a factor e~^, as de- 
scribed in equation (|4|. Decreasing Zd therefore increases the den- 
sity gradient, and the density ratio between the east and west sides 
of SS 433. Increasing the density gradient means that SS 433's west 
jet decelerates rapidly due to the sharp increase in density on the 
west side of SS 433, in contrast to the east jet, which decelerates 
much less due to a rapid decrease in the ISM density over short 
distances to the east of SS 433. 

A series of jet simulations were performed to test many com- 
binations of {no,Z^), in order to reproduce the correct east- west 
asymmetry ratio observed in W 50. These results are presented in 
Figure |6] from which two parameter combinations appear to repro- 
duce W50's east- west lobe extents very well, as shown in Table |4] 



7.2 Collimation of the jet lobes 

It is important to realise that the two best-fit parameter sets in|4]are 
only relevant for a jet mass-loss rate of Mjet = lO'^^Mgyr'^ This 
realisation comes from a particularly important result demonstrated 
by the dashed red- white and orange- white lines in Figure [6] corre- 
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Table 4. Best-fit parameter values. 



Parameter Set 


Wo (particles cm ^) 


Zd (pc) 


Setl 


0.2 


40 


Set 2 


0.1 


30 



spending to simulations Lobeslb and LobesSb respectively. These 
simulations show that changes in the background density normali- 
sation riQ and the jet mass loss rate Mjet can become indistinguish- 
able, as previously hypothesised in ^ For example, increasing the 
jet power by a factor of 10 (equivalent to increasing Mjet if the jet 
speed is a constant) whilst keeping all other parameters constant, is 
equivalent to reducing the background density normalisation (the 
resistance felt by the jets) by a factor of 10. Thus, I define the quan- 
tity to indicate the penetrative ability of a given jet: 



such that any two jets with the same value of in will follow 
the same path in Figure |6] for a given value of Zd. Thus, for any 
given mass-loss rate, the corresponding best-fit value for uq can be 
calculated from ( [48] ). 

7.3 Hydrodynamic refocusing of SS 433 's jets. 

It has been suggested that astrophysical jets of the hollow conical 
type, may undergo refocusing towards the centre of the cone ( |Eich-| 
|ler|1983 ) due to the ambient pressure from the ISM through which 
the jet propagates. However, we find evidence of refocusing of our 
hollow conical jets through a diff'erent mechanism, which is purely 
dependent upon the kinematics of the ISM within the cocoon of 
the jet. Two subtly diff'erent refocusing mechanisms are suggested 
here, but both are based upon momentum exchange between the jet 
and its environment, as follows: 

(i) Kinetic refocusing - the jets are hydrodynamically refocused 
via collisions with the low-density but fast-moving and turbulent 
ISM within the jet cocoon, which has been energised through past 
interactions with earlier jet ejecta. Due to the low density of the gas 
within the jet cavity, it is quite easy for this material to be acceler- 
ated to high speeds through interactions with the jets. After gaining 
momentum from the jets, the fast-moving gas interacts with the 
wall of the cocoon and develops a vortex-like rotation, as indicated 
by the arrows of the velocity field of the gas superimposed onto Fig- 
ure[7^. In this case, momentum transfer from the jets to the cocoon 
ISM has the dual efi'ect of altering the jet trajectory and further stir- 
ring the turbulent ISM motion. This eff'ect is self-sustaining once it 
begins. This eff'ect is seen in all of our simulations that have a con- 
stant jet cone angle greater than zero, such as Jet 5, Jet 6 and Jet 7 
of Figures [5^, [Sf and[5^. The process of momentum exchange with 
the turbulent ISM is depicted in Figure [7]3. 

(ii) Static refocusing - a more gradual focusing mechanism oc- 
curs when the jets interact at a grazing angle with the dense, but 
comparatively slow-moving, cavity walls of the jet cocoon. The 
interaction occurs at a grazing angle when the jet cocoon is suf- 
ficiently elongated, such as that created by a cylindrical jet (or a 
jet with a near-zero cone angle). Due to the high density of the 
cavity wall, it experiences a much smaller acceleration per unit 
momentum-exchange than does the low-density gas of the kine- 
matic refocusing method. As a result, the jet is redirected back to- 




log[/9 (gcm"^) ] -25.3 -25.0 -24.7 -24.4 -24.1 -23.8 -23.5 -23.2 



Figure 7. Demonstration of the kinetic hydrodynamical refocusing of hol- 
low conical jets. The small black velocity vectors in (a) have magni- 
tude proportional to the square-root of the average speed at each point, 
and are overlaid upon a zoomed-in density map of Jet 7. The large 
white arrows in (a) show the net vortex-like movement of the turbu- 
lent gas in each quadrant. These vortices tend to move around appre- 
ciably within each quadrant, which is the reason for the asymmetry in 
their locations in this snapshot. The black arrows in (b) indicate the ini- 
tial trajectory of the jets. The purple arrows indicate the turbulent mo- 
tion of the agitated gas within the jet cocoon, and the red arrows in- 
dicate the subsequent interaction and refocusing of the jets back to- 
wards the jet-axis. A short animation showing this effect is avilable at: 
\http : / /www-astro. phy sics.ox.ac.uk/ ~ptg/RESEARCH/research.html\ 

wards the mean jet axis upon collision with the cavity walls. This 
mechanism is seen to occur in Jet 8 of FigureBJi. 



8 SNR + JETS: RESULTS AND DISCUSSION 

Following from Table |4] in ^7.1| the combination of parameters 
(no - 0.2cm~^,Zd - 40 pc) and (hq = 0.1 cm~^,Zd - 30 pc) were 
determined to be most representative of the observations of W 50, 
for a jet mass loss rate of Mjet = 10~^Moyr~^ Although both of 
these parameter combinations are equally feasible, they are also al- 
most equivalent in terms of the jet morphology produced. Hence, 
we focus on the first parameter set (hq = 0.2cm~^,Zd = 40 pc) as 
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Table 5. A description of the key parameters relevant to each of the SNRJet simulations. 



Shortname 


Description 


Parameter Details 




SNRA 


SNR only 


SNR defaults and: p-profile = 1, v-prof =1, no = 0.1 cm~^, Zd 


= 40pc 


SNRB 


SNR only 


SNR defaults and: p-profile = 1, v-prof = I, no = 0.2 cm~^, Zd 


= 30pc 


SNRJetl 


SNRA with a modified Jetl 


Modified parameters: no = 0.1 cm~^ 




SNRJet2 


SNRB with a modified Jetl 


Modified parameters: Zd = 30 pc 




SNRJetS 


SNRA with a modified Jet 6 


Modified parameters: = 20° 




SNRJet4 


SNRA with a modified Jet 8 


Modified parameters: 6* = 20°/1500jr 




SNRJet5 


SNRA with a modified Jet 9 


Modified parameters: 01=0°, 02 = 20° 





Jet default parameters: Mjet = 10 ^ M© yr ^ 



being the best description of SS 433 's environment, and the effects 
of the interaction between the jets and SNR are investigated for this 
setting. 

All jet models featured here use the same jet speed of 0.26c 
corresponding to the observed speed for SS433's current jet state. 
We acknowledge that the jet speed is one of the variable parameters 
that can change between outbursts of jet activity, and the choice to 
keep the jet speed constant here has no physical basis other than 
to limit the number of free parameters in the models. This enables 
us to focus on the effects of changing the jet precession cone angle 
with each different jet outburst. The effects of varying other param- 
eters such as jet speed and jet mass loss rate, will be investigated in 
a subsequent publication. 



8.1 Simulating the Jet-SNR interaction 

The best-fit parameter sets from Table [4] were used to create two 
new SNR simulations: SNRA and SNRB, thereby incorporating the 
most appropriate background environment settings for the evolu- 
tion of the new SNR shells. These new SNR simulations are oth- 
erwise identical to SNR 6. As per j4.1| the SNR evolution is fol- 
lowed until the radius approaches 45 parsec, which happens at 
^21,000yrs for SNRA, and ^17,000yrs for SNRB, due to the 
slightly lower density of SNRB. Various jet models were then in- 
voked to investigate the SNR-jet interaction, and parameters tested 
in each simulation are described in Table [5] 



8.2 Describing each SNR- Jet interaction 

The final snapshot for each of the simulations described in Table 
|5]are shown in Figure [8] where the left column shows the density 
maps, and the right column shows the gas speed (v^ + v^)2. The 
purple contours on the left column of Figu re [8] show the outline 
of W 50 from the VLA radio observations of |Dubner et al.| ( p^98] ). 
The green contours on the right column of Figure[8]show the struc- 
ture of SS433's jet lobes from the ROSAT X-ray observations of 
[Brin kmann e t al.| ( [l996| ). The inlay on each panel shows a 4x mag- 
nified region about the coordinates of the jet launch. 



the purple contours in Figure [Sj a). This confirms that the parame- 
ters were determined correctly (no = 0.2 cm~^, Zd = 40 pc), and that 
the presence of the SNR has little effect upon the east-west lobe 
asymmetry along the jet-axis direction. The SNR does however 
widen the jet at the point where the jet penetrates the SNR shell, 
in a similar manner to waves diffracting through an aperture. This 
is because the jet breaks through the relatively high density of the 
compressed SNR shell into the surrounding lower-density ISM, and 
the impact of the jets with the shell increases the momentum of 
the SNRJet material in the direction perpendicular to the jet axis. 
Note that the jet lobes produced in SNRJetl have a higher degree 
of collimation than those of Jetl from Figure [5^, due to the back- 
ground density normalisation being a factor of 5 lower in SNRJet 1 
than Jet 1 . Judging by the purple contours, the east lobe of SNRJet 1 
seems to be approximately equal in width to the eastern radio lobe, 
whereas the transition between the SNR shell and the western lobe 
seems to be more gradual in the purple radio contours than the sim- 
ulation would indicate. This could be due to a previous jet episode 
impacting on the younger SNR at an earlier time, such that the SNR 
shell is inflated more rapidly at the location of the jet impact. The 
green contours in FigurejSjb) indicate again that the eastern jet lobe 
has similar dimensions for both the simulation and X-ray observa- 
tions, however the western lobe X-ray emission is abruptly trun- 
cated short of the full extent of the simulated jet lobe. The X-ray 
lobe is both wider and shorter than the west lobe of the simulation, 
as though the X-ray emitting jet material is hitting a dense wall of 
material, which is true if the increased radio brightness of the neb- 
ula at this position indicates an increase in density (see point D in 
Figure pT|). Again, this could be indicative of a previous jet ejection 
episode. 

SNRJet 2 of Figure [Sjc & d) is similar to SNRJetl in most 
respects, except for the fact that the jet lobes are more highly colli- 
mated than SNRJet 1 due to a decrease in the density normalisation 
by a factor of 2. The jet lobes produced by SNRJet 2 are thinner 
than the radio/X-ray lobes as indicated by the purple/green con- 
tours, and this probably means that the density of no = 0.1 particles 
cm~^ is too low, and it's likely that SNRJet 1 with no = 0.2 particles 
cm~^ is more representative of the environment at the locality of 
SS433. 



8.3 The roles of ambient density no and scale-height Zd 

SNRJet 1 of FigurejSja & b) shows how Jet 1 from Figure[5^ would 
interact with SNRA. The extents of the east and west jet lobes are in 
good agreement with the radio observations of W 50 as shown by 



8.4 The role of precession cone opening angle 

SNRJet 3 of FigurejSje & f) features a precessing jet with half-cone 
angle of 20°, very similar to that of the current jet in SS433. Al- 
though this jet model doesn't reproduce the lobes of W50, it does 
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Figure 8. This figure illustrates the varying nebular morphologies produces through the interaction of 5 different jet models with an evolved SNR. Refer to 
Table [5] for a summary of the differences between these simulations, and ^for a detailed commentary on the interpretation of each one. The animations 
associated with these simulations are available at: http : / /www-astro. physics. ox.ac.uk/~ptg /RES EARCH/ research.html 
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Figure 9. The image shows a density map (logarithmic colour scale, 
white = high density, black = low density) of SNR A in the final stages of 
evolution. The blue line shows a density line profile taken through the centre 
of the displayed image (indicated by the dotted fine). The green fine shows 
the same density line profile of the unperturbed ISM (before the supernova 
explosion occurred). This figure illustrates the diff'erence in the local den- 
sity distribution as seen by Jet 8 (green line, without SNR), and SNRJet4 
(blue fine, after SNR has ploughed through the ISM). The lower density 
environment through which the jet from SNRJet4 propagates, is the reason 
for the morphological diff'erences with Jet 8. 



accelerate the expansion of the east and west hemispheres of the 
SNR bubble. This is further evidence that a previous jet episode 
with a large cone angle probably contributed to the shape of W 50. 
It is unlikely that this inflation of the east and west hemispheres 
were caused by the current jet episode in SS 433 due to time argu- 
ments, as explained in ^8.6| This hollow conical jet model is also 
affected by the Kinetic Refocusing mechanism described in g7.3| 
and the effects are noticeable at around 500 years onwards in the 
simulation (refer to the simulation movies at the accompanying we- 
blink). 



8.5 The role of time-varying precession cone opening angle 

SNRJet4 of Figure [8|g & h) features a conical jet whereby the jet 
precession cone half-angle linearly increases from 0° to 20° over a 
period of 1500 years. The jet parameters are almost identical to that 
of Jet 8 from ^ with the only difference being the time taken for 
the jet to reach the maximum cone angle of ~20°, which was 1000 
years rather than 1500 years. However, the difference in the nebular 
morphology produced by Jet 8 and SNRJet4 is striking (compare 
Figure [5ji with Figure [8^), and this difference is due entirely to the 
presence of the compressed SNR shell. 

The hypersonic Shockwave created by the supernova explo- 
sion sweeps up most of the gas in the ISM, and the density distri- 
bution from the final snapshot of SNR A is shown in Figure|9] along 
with the original density profile of the background ISM before the 
supernova explosion. The morphological difference between Jet 8 
and SNR Jet 4 demonstrates (yet again) the importance of the ISM 
density distribution in the vicinity of the jet launch site, and the 
profound influence this can have upon the resulting jet cocoon. 



8.6 The role of episodic jet ejection 

SNR Jet 5 of FigurejSji & j) uses an episodic jet model with two out- 
bursts; the first jet outburst features a cylindrical jet which lasts for 
200 years, and the second outburst is a precessing conical jet with 
half-angle 20° of 2300 years in duration. SNRJetS displays both 



Kinetic and Static refocusing (see the movies on the accompany- 
ing webpage). Kinetic refocusing becomes noticeable around ~700 
years after the jet is first switched on, and is readily recognisable 
as four regions of swirling gas within the SNR, which are (approx- 
imately) symmetrically located about the jet launch point (one in 
each quadrant). Static refocusing begins approximately 2000 years 
after the jets are switched on, whereby the jets are refocused to- 
wards the jet axis upon interaction of the denser jet cocoon gas 
located near the SNR shell, which was penetrated several hundred 
years earlier by the cylindrical jet. 

This jet model was devised to combine the successes of both 
cylindrical jets SNRJetl in reproducing the correct jet lobe extents, 
and the conical jets SNRJetS in producing the required inflation of 
the east and west sides of the SNR such that a gradual transition 
from SNR shell to lobes is created. However this attempt was un- 
successful, because the cylindrical jets (despite the short episode 
duration of just 200 years) exceed the extent of W 50's lobes in the 
time it takes for the conical jets to catch up with the SNR shell and 
contribute to its expansion, as shown in Figure [Sj. 



8.7 Jet Kinematics 

Astrophysical jets sweep up the ISM just as supernova blast waves 
do, with the very important difference that supernova explosions 
are very short-lived (in astrophysical terms), injecting enormous 
amounts of energy into the ISM in one event. For jets, the ejection 
of high-speed material often continues for long periods of time, and 
in fact many thousands of years of jet activity are required to create 
the W50 nebula surrounding SS 433. As a consequence, the prop- 
erties of the environment surrounding the jet-launch site, change 
with time. The activity of the jets clears a path through the ISM 
for the jet material that is ejected at later times. This is one pos- 
sible explanation for the variable observational characteristics of 
microquasar jets seen at different times or outbursts, as mentioned 
in g2.3| However, a more obvious result from these simulations is 
that the jet suffers very little deceleration inside the jet cocoon or 
cavity produced by earlier jet activity, and most of the decelera- 
tion occurs in the final stages when the jet ejecta catches up with 
the dense build-up of gas behind the jet shock front. This is illus- 
trated in the right-hand column of Figure [8] where the colour maps 
show the gas speed; the pink/white regions show gas travelling very 
close to the launch-speed. The full details of this effect are given in 
a companion paper (how to site this companion paper?). 



8.8 Jet symmetry and collimation 

The level of collimation is actually higher for pulsed jets than for 
continuous jets. For a continuous jet, although the jets receive con- 
stant replenishment of the kinetic energy and therefore don't show 
signs of deceleration, a side-effect is that the jet bunches-up and 
crashes into the jet in front. Any tiny asymmetric displacement of 
the bunched-up jet about the jet axis will cause the upstream jet 
to transfer an increased amount of momentum in the transverse di- 
rection. An analogy for this effect is drawn with that of a cue-ball 
struck off-centre by the cue. This effect is much less noticeable 
when the jet is pulsed, presumably due to a lower occurrence of jet 
blockages. 



© 0000 RAS, MNRAS 000, 000-000 



20 P, T. Goodall, E Alouani-Bibi and K. Blundell 




I— ' ' ' " ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' >- 

-120 -90 -60 -30 30 60 90 

Distance from SS 433 (pc) 



Figure 10. The vertical width of each nebula produced from the SNR-jet 
interaction simulations from Figure [8] is shown as a function of distance 
from SS433 by the coloured lines in each panel. The width of W50 as a 
function of distance form SS 433 is also shown by the thick white line in 
each panel, as measured from the Dub ner etal.H1998) image. Note how the 
transition between SNR and jet lobes is very smooth in reality (white lines), 
and in fact it is difficult to say where the SNR shell ends and the jet lobes 
begin when viewed in this way. The dashed coloured lines show the width 
of the jet internal to the SNR. The bumps in the lobes indicate the locations 
of the annular features in W 50's lobes. 



8.9 Limitations of current physical models 

Due to the microphysics not being well understood, comparison be- 
tween observation and simulations are not straightforward. Since 
we are concerned with the general morphology of the W 50 neb- 
ula, only the overall brightness profile from the radio observations 
should be compared with Fig ure [5] and Figure [8] Optical observa- 
tions of W50 ( |Boumis et al.||2007^ ) show small-scale filamentary 
structures; simulations of these filaments would require a detailed 
description of the relevant microphysics (such as the amount of 
electron heating in collisionless shocks, the detailed structure of 
SN ejecta). We note that the presence of any clumpyness of the 
circumstellar medium and possibly dust would have an eff'ect on 
the propagation of the young SNR and jets. The density maps from 
our simulations reproduce the key morphological aspects of W50 
such as the lobes with appropriate lengths and widths, and a cir- 
cular region away from the jets due to the expansion of the SNR, 
in good agreement with the radio observations. The X-ray obser- 
vations appear to trace out regions of recent jet activity, but again 
only the general large-scale shape is important for comparison with 
our simulations. 



9 CONCLUSIONS 

The cylindrical jet models reproduce the lobes of W50 very well 
in terms of their absolute and relative extents in the east and west 
directions. As Figure[8]and Figure [T0| show, the cylindrical jets cre- 
ate jet lobes with well-defined boundaries at the site of the impact 
between jets and the SNR shell. This is not in agreement with the 
radio observations of [Dubner et al.| ( p^98] ) in which there is a more 
gradual transition between the circularity of the SNR shell and the 
elongated lobes (see Figure [TO]). 

The conical jet models do produce lobes with a lower degree 
of collimation appropriate to W50's lobes very near to the SNR 
shell, however these jets continue to diverge away from the jet axis 
and prove too wide at larger distances from SS 433 to be considered 
an accurate representation of W 50's lobes. 

It was hypothesised that the smooth transition between SNR 
shell and lobes might be created if a conical jet model is used in 
conjunction with a cylindrical jet model. This was tested in SNR- 
Jet5 where a cylindrical jet (an earlier jet episode) was invoked 
before a conical jet with the precession cone half-angle approxi- 
mately equal to that currently observed in SS433. However, this 
leads to problems due to the large time interval (at least 2000 years) 
required for the conical jets to reach the surface of the SNR shell, 
during which time the cylindrical jet travels beyond the boundaries 
of the simulation domain, indicating that a real jet would have far 
exceeded the extent of W 50's lobes in the time required by the con- 
ical jets. This problem is eradicated if the conical jets begin first and 
are allowed to evolve for as long as is necessary before the cylindri- 
cal jets begin, because the latter only require of order ~1600 years 
to reach the required lobe extent. However the situation gets more 
complex in that the final jet episode must be in agreement with cur- 
rent observations, and so a third jet episode representing SS433's 
current jet precession state must be invoked in order to satisfy this 
requirement. The role of episodic jet outbursts will be further in- 
vestigated in future work. 

In summary, it seems clear that neither cylindrical nor conical 
jets can independently reproduce the interesting morphology dis- 
played by W 50, and it seems several jet episodes from SS 433 with 
varying jet characteristics are needed in order to sculpt W 50 in the 
ISM. 



9.1 Energetics of the SS 433-W 50 system 

With a radius of 45 pc, W 50 is among the largest SNR observed 
to date, and to consider higher SNR blast energies of 10^^ ergs or 
so, would be reasonable. This would usefully lessen the latency 
period between the supernova event and the jets switching on, since 
the SNR expansion speed begins as ^JlE^^^^JM^y Such high blast 
energies would be expected for a very massive progenitor. 

If both stars in the binary formed at the same time, a simple 
main-sequence lifetime argument requires that the SS433's pro- 
genitor must have been more massive than its Wolf-Rayet compan- 
ion star at 24Mo |Blundell et al.t 2008 ), since it was first to detonate. 
A massive Wolf-Rayet progenitor producing a Hypernova explo- 
sion includes the possibility that SS 433 may have formed through 
a Gamma Ray Burst (GRB) event. To simulate this scenario, our 
models would need to be adapted to include a short but highly rela- 
tivistic jet outburst simultaneous with the progenitor detonation, to 
be investigated in a further publication. 

The rapidly expanding SNR could make W 50 an ideal candi- 
date for producing of high energy cosmic rays. We note that with 
the large radius of 45 pc for W50 and the Galactic field strength 
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Figure 11. The striking similarities and differences between (a) the VLA 
radi o image of W50 (Dubner et al. J998; > and (b) the ROSAT X-ray im- 
age jBrinkmann et al.|l996) are entirely consistent with the hypothesis of 
episodic jet ejection in SS 433. The ROSAT contours are overlaid upon the 
radio image in (a), and the radio contours are overlaid upon the ROSAT 
image in (b). 



of 5yuGauss makes confinement of particles with energies up to 
2xl0^^eV possible. 



9.2 The Radio and X-ray correlation 

Our simulation results can be shown to be consistent with radio and 
X-ray observations of W50 as shown in Figure [TT] in conjunction 
with the jet outbursts hypothesis. Note the striking correspondence 
in Figure [TTJb) between the low-level (blue) emission and the out- 
line of the contour from the radio nebula. Furthermore, the bright 
filament at point A on the far west of W 50 is coincident in both 
X-rays and radio. It is equally interesting to note that the brightest 
parts of the X-ray lobes (emission from the jets) at points B and C 
are almost exactly coincident with regions where the radio emis- 
sion is faintest (i.e. radio-quiet regions where the X-ray lobes are 
brightest). Finally, the point labeled D shows an abrupt termination 
of the eastern X-ray jet lobe followed by an abrupt brightening of 
the radio jet lobe. This could indicate two diff'erent populations of 
particles, perhaps where the relatively recent, hot, and X-ray emit- 
ting jet ejecta from a new jet episode has caught up with the older, 
cooler, ejecta from an outburst that happened much further in the 
past. 

It is generally noted that the X-ray emission from W 50 is seen 
to be confined to a region much closer to the mean jet axis than the 
precession cone of SS 433 's jets would indicate, and this could be 
attributed either to a previous jet ejection episode with a preces- 
sion cone angle close to zero, or to the mechanism of refocusing of 
SS 433's jets as detailed in fJ3\ 



9.3 Comparison with previous work 

Despite the appeal of some of the previous hydrodynamic mod- 
els attempting to reproduce W50, there are some inconsistencies 
which must be discussed in order to begin to fully understand the 
nature of the SNR-jet interaction in SS433 and W50. First, we 
note that concrete graphical evidence in the literature supporting 
claims of having reproduced the annular structure in the lobes of 
W50 through hydrodynamic simulations, is (so far) scant. In Fig. 
3 of Velazquez & Raga (2000 ), the pressure contour at 1150 years 
shows the emergence of two annuli through the jet-SNR interac- 
tion, however it is not clear if these are artefacts of the reflective 
boundary conditions imposed. The computed synchrotron emission 
(Fig. 5 of Velazquez & Raga (2000)) does resemble the annular 
structure in the eastern lobe of W 50 at approximately 3000 years 
after the jets are initiated. However, it is easy to create artificial 
annular structures through azimuthal revolutions with axial sym- 
metry. The follow-up paper |Zavala et al.| |2008) with improved jet 
models and in three dimensions, reports no evidence of any annu- 
lar structure. They acknowledge that the annular structures created 
in the previous paper are a possible consequence of the reflective 
boundary conditions imposed therein, and the undesirable eff'ects 
produced via axisymmetric codes with reflective boundary condi- 
tions are well known to those authors { [Raga et al.|2007| ). 

Each of the simulations shown so far in the literature have 
featured scaled-down versions of the W50 system. The primary 
assumption for any work involving the SS 433-W 50 system, is that 
both objects are at the same location, and independent methods 
have determined a distance of 5.5 kpc is appropriate to each of 
SS 433 dBlundell & Bowler|2004 ') and W 50 (Lo ckman et al.|2007| l. 
At this distance, the radius of the SNR is 45 parsecs and the jet ex- 
tents are 86.5 and 121.5 parsecs for the west and east jet lobes re- 
spectively, as indicated in Figure [T] These physical sizes (and their 
ratios) are of paramount importance in forming a consistent model 
of the SS 433-W 50 interaction, that is in agreement with obser- 
vations. We stress that (in contrast to the claims of [Zavala et al.| 
([2008)), it is not sufficient to simply scale-up simulated models of 
the system for comparison with observations. The reasons for this 
are as follows: 

(i) The absolute extents of the jet lobes at a given time are char- 
acteristic of the integrated column density along the path of the jets 
(a measure of the resistance faced by the jets) as they travel through 
the ISM. Similarly, the radius of the SNR at a given time is charac- 
teristic of the mass it has swept up, which depends roughly on the 
cube of the radius. Aside from the density normalisation, the evo- 
lution of each of these is also dependent upon the energy input, as 
governed by ^biast for the supernova, and Mjet over some time pe- 
riod for the jets. The jets propagate much faster than the SNR at all 
times, and so a morphological match to W 50 at some scaled-down 
radius of the SNR will be lost at later times. In essence, the jets 
would extend far beyond the domain of W 50 by the time the SNR 
has expanded to a radius that is in agreement with observations. 

(ii) The relative extents of the east and west jet lobes are charac- 
teristic of the exponential density background, and the ratio of the 
lobe extents is not linearly scalable with either total jet travel time 
or total lobe extent, as shown in Figure |6] 

(iii) Simulating the absolute extents of both SNR and jet lobes 
is important as this allows us to place constraints on the formation 
time and energetics for each component, and has relevance to the 
aforementioned latency period. 

(iv) Jets could possibly inflate (accelerate the shell expansion of) 
a small SNR corresponding to some scaled-down version of W50. 
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However, we find that this is not the case for a SNR as large as 
45 pc, and the jets have absolutely no eff'ect upon the SNR expan- 
sion in the direction perpendicular to jet launch. 

9.4 Summary 

Of all simulations presented in this paper, only cylindrical jets have 
produced sufficiently collimated jet lobes comparable to the jet 
lobes in W50. Simulations mimicking the jet precession cone ob- 
served in SS 433 at the current epoch, do not reproduce the lobes of 
W 50. Neither cylindrical, nor conical jets resembling SS 433's cur- 
rent precession state, are able to reproduce the interconnecting re- 
gion between W 50's circular SNR shell and the jet lobes, although 
very large cone angles of 40° can produce it. Our findings there- 
fore suggest that at least three very different jet states are required 
during SS 433's very active jet history, in order to produce (i) well 
collimated jet lobes, (ii) smoothly connects SNR and lobes, (iii) the 
precessing corkscrew jet observed at the present day. This is the 
minimum number of jet outbursts required to create a W 50 shaped 
nebula around the microquasar, but it is possible than many more 
jet outbursts have made contributions over time. Further evidence 
of episodic jet outbursts from SS 433 comes from the comparison 
of the radio and X-ray observations detailed in |9.2| and from an 
archival radio study of W 50's kinematics described in a companion 
paper (cite the companion paper), the results of which are entirely 
consistent with the hydrodynamic simulations presented here. 

The resultant nebula morphology created through sequential 
episodic jet activity is strongly dependant upon the nature of the 
jet activity that occurred previously to it. Early jet activity ploughs 
gas away from the point of jet-launch, thereby creating a cavity 
of lower ISM density through which later jet ejecta can travel more 
easily. Nearby ISM that has been agitated or stirred by the influence 
of the jets, can create vortices and bulk flows of gas which kineti- 
cally refocus conical jets, through exchange of momentum. This jet 
refocusing also occurs when the jets interact with the dense cavity 
walls behind the jet-lobe shock front. 

The focal length associated with the kinematic refocusing 
mechanism of ^7.3| increases slightly with the precession cone an- 
gle and the size of the jet cavity, but is typically of order 10 pc or so. 
The focal length associated with the static refocusing mechanism of 
^7.3| also depends upon the precession cone angle, and the width of 
the jet cavity carved out by the earlier jet ejecta, but at 30-40 pc 
it is typically much larger than the focusing length induced by the 
kinematic mechanism. As such, the focusing length for SS 433's 
jets should tell us something about the history of its jet activity, 
and an investigation is under-way to observe this. We stress how- 
ever, that while these are both interesting hydrodynamical efi'ects, 
neither method of jet refocusing produces the level of collimation 
required to describe the lobes of W50, based on SS 433's current 
jet precession state. 

We further conclude that careful hydrodynamical simulations 
can provide useful constraints upon the parameters of complex as- 
trophysical systems. The distinct morphology of the W 50 complex 
has allowed us so determine the following parameters as being most 
appropriate in describing SS 433's environment: {uq = 0.2cm~^, 
z - 40pc) and (no - O.lcm"^, z - 30pc). 
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